Thermal vortex dynamics in thin circular ferromagnetic nanodisks 
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The dynamics of gyrotropic vortex motion in a thin circular nanodisk of soft ferromagnetic mate- 
rial is considered. The demagnetization field is calculated using two-dimensional Green's functions 
for the thin film problem and fast Fourier transforms. At zero temperature, the dynamics of the 
Landau-Lifshitz-Gilbert equation is simulated using fourth order Runge-Kutta integration. Pure 
vortex initial conditions at a desired position are obtained with a Lagrange multipliers constraint. 
These methods give accurate estimates of the vortex restoring force constant kp and gyrotropic 
frequency, showing that the vortex core motion is described by the Thiele equation to very high 
precision. At finite temperature, the second order Heun algorithm is applied to the Langevin dy- 
namical equation with thermal noise and damping. A spontaneous gyrotropic motion takes place 
without the application of an external magnetic field, driven only by thermal fiuctuations. The 
statistics of the vortex radial position and rotational velocity are described with Boltzmann distri- 
butions determined by kp and by a vortex gyrotropic mass mc ~ /kp, respectively, where G is 
the vortex gyrovector. 

PACS numbers: 75.75.-c, 85. 70. Ay, 75.10.Hk, 75.40.Mg 
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INTRODUCTION: VORTEX STATES IN 
THIN NANOPARTICLES 



Vortices in nanometer-sized thin magnetic particlesi 
have attracted a lot of attention, due to the possibili- 
ties for application in resonators or oscillators, in detec- 
tors, as objects for data storage?? We consider the dy- 
namic motion of an individual vortex in a thin circular 
disk (radius R and height i <C i?) of soft ferromagnetic 
material such as Permalloy-79 (Py) where vortices have 
been commonly studied.— In disks of appropriate size, 
the single vortex state is very stable and of lower en- 
ergy than a single-domain statci^ Especially, we study 
the effective force constant kp responsible for the restor- 
ing force on a vortex when it is displaced from the disk 
center, F = — fci?X, where X is the vortex core position 
relative to the disk center. If the vortex is initially dis- 
placed from the disk center, say, by a pulsed magnetic 
field, ^ it oscillates in the gyrotropic mode,^ at an angu- 
lar frequency log = kp/G, where G = Gz is the vortex 
gyrovector, pointing perpendicular to the plane of the 
disk. The vortex gyrotropic motion has been observed, 
for example, by photo emission electron microscopy using 
x-rays 4 Not only in small disks but in many easy-plane 
magnetic models this type of vortex dynamics has been 
studied for its interesting gyrotropic dynamicsi^iiS The 
gyrotropic mode is due to the translational modeii*^ in 
a whole spectrum of internal vibrations of a magnetic 
vortexJ^ 

The force constant estimated analytically by Guslienko 
et al^ using the two-vortices mode]~ gave predictions 
of the gyrotropic frequencies obtained in micromagnetic 
simulations, with reasonably good accord between the 



two. Here, we discuss direct numerical calculations of 
kp based on static vortex energies, by using a Lagrange 
multipliers technique to secure the vortex position X at 
a desired location, and map out its potential within 
the disk. We apply an adapted two-dimensional (2D) 
micromagnetics approach for thin systems to calculate 
the demagnetization field. As a result, the calculations 
can be directly compared with the two- vortices prediction 
for kp. 

At the same time, we give a corresponding study of 
the vortex dynamics to calculate the gyrotropic frequen- 
cies. It is found that the static results for kp combined 
with the dynamics results for log agree with the predic- 
tion of the Thiele dynamical equation^i^iii wg = kp/G, 
to very high precision. Similar to Ref. 0, kp is found 
to be close to linear in the aspect ratio L/R if the disks 
are thin but not so thin that the vortex would be desta- 
bilized. Our values for kp are slightly less than those 
in the two-vortices model, as the numerical relaxation of 
the vortex structure allows for more flexibility than an 
analytic expression. 

A micromagnetics study of this systemic for finite tem- 
perature shows evidence for a spontaneous gyrotropic vor- 
tex motion with a radius of a couple of nanometers, with- 
out the application of a magnetic field. The spontaneous 
gyrotropic motion occurs even if the vortex is initiated 
at the center of the disk in the simulations. It is clear 
that thermal fluctuations should lead to a random dis- 
placement of the vortex core away from the disk center, 
however, it is striking that the ordered gyrotropic rota- 
tion appears and even dominates over the thermal fluc- 
tuations. Here we conflrm this effect, and also find that 
a spin wave doublelfi^ (of azimuthal quantum numbers 
TO = ±1) is excited together with the gyrotropic motion. 
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Having at hand the force constant kp, we can analyze 
both the dynamics and the statistics of the gyrotropic 
motion induced by the temperature. The study of the 
finite temperature dynamics is carried out using a mag- 
netic Langevin equation that includes stochastic mag- 
netic fields together with damping. We discuss the so- 
lution via the second order Heun method^°'^^ applied to 
magnetic systems. Further, we introduce a technique for 
estimating the location of the vortex core accurately in 
the presence of fluctuations. Based on the behavior of 
kp with disk geometry, we find it possible to predict the 
RMS displacement of the vortex core in equilibrium. By 
using the collective coordinate Hamiltonian for the vor- 
tex, as derived from the Thiele equation, it is also pos- 
sible to determine the probability distributions for vor- 
tex radial displacement r = |X| and rotational velocity 
V = ujgt- It is interesting to see that the velocity dis- 
tribution, f{V), is of the Boltzmann form for a particle 
with an effective mass given by mc — G"^ /kp, which is 
found to depend only on the gyromagnetic ratio 7, the 
magnetic permeability of free space /iq and the disk ra- 
dius. 



II. DISCRETE MODEL FOR THE CONTINUUM 
MAGNET 

We determine the magnetic dynamics for a continuum 
magnetic particle, but using a thin-film micromagnet- 
ics approach, ^=2. defining appropriate dipoles at cells of 
a two-dimensional grid. This is a modification of usual 
micromagnetics^- where a 3D grid is used. The parti- 
cle has a thickness L along the z-axis, and a circular 
cross-section of radius R. For thin film magnets it is rea- 
sonable to make the assumption that the magnetization 
M(r) does not depend on the coordinate z through the 
thickness. This is acceptable as long as the particle is 
very thin. The demagnetization field tends to cause M 
to lie within the xy plane in most of the sample^ except 
for the vortex core region. Even in the vortex core, how- 
ever, one should not expect large variations of M with z, 
due to the dominance of the ferromagnetic exchange over 
the dipolar interactions through short distances. In this 
situation for very thin magnets, this 2D approach has the 
obvious advantage of greater speed over 3D approaches, 
without sacrificing accuracy. It is somewhat like using a 
single layer of computation cells in 3D micromagnetics, 
with the cell height longer than its transverse dimensions. 

The energy of the original continuum system, including 
exchange and magnetic field energy, can be expressed as 
a volume integral. 
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The magnetization scaled by saturation magnetization 
Ms is used to define the scaled magnetization m = 
M/Ms, that enters in the exchange term, where A is 



the exchange stiffness (about 13 pJ/m for Permalloy). 
The last term is the interaction with an externally gen- 
erated field, H°'^*^. The demagnetization energy involves 
the demagnetization field that is generated by M, 
and which is determined through a Poisson equation in- 
volving the scalar magnetic potential 4>m, 
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(2) 



This is solved formally in three dimensions using a con- 
volution with the 3D Green's function: 

$M(r) - I dh' GMr-r')pM{r'), (3) 
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(4) 



However, this is reduced to an effective 2D Green's op- 
erator, appropriate for thin magnetic film problems, re- 
viewed below. 



A. The micromagnetics model 

The micromagnetics^^i^ is set up to use nanometer- 
scaled cells in which to define coordinates as the av- 
eraged scaled magnetization in that cell. The system is 
divided into cells of size a x a x L (axais the cross 
section in the xy-plane), rather than cubical cells. Each 
cell i contains a magnetic moment of fixed magnitude 
fi — La? Ms, where Ms is the saturation magnetization. 
The direction of the (assumed uniform) magnetization 
in a cell is a unit vector, rhi, whose dynamics is to be 
found. The cells interact with neighboring cells via the 
exchange interaction, and with all other cells, due to the 
demagnetization field, and also with any external field. 

For the square grid of cells, the exchange energy is 
found to be equivalent to 
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(5) 



where the sum is over nearest neighbor cell pairs. The 
energy scale of exchange is taken as the basic energy unit. 
Thus it is convenient to define an effective exchange con- 
stant acting between the cells, 

J = 2AL, (6) 

and for the computations, all other energies will be mea- 
sured in this unit. In addition, the saturation magneti- 
zation is a convenient unit for magnetic fields as well as 
for M. So we define scaled fields, 
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(7) 



As a result of this, the magnetic field interaction energy 
terms are scaled here as follows. For the demagnetiza- 
tion. 
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and for the energy in the external field, 



mi 



(9) 



These depend on the definition of the exchange length, 
\o. = ^l^^, (10) 



/ 2A 



that gives a measure of the competition between ex- 
change and dipolar forces. This means that the effective 
2D Hamiltonian can be written as 



H 




(11) 



B. The demagnetization field H'^' in a thin film 

It is important to calculate the demagnetization field 
efficiently and accurately, as it plays an important role 
in the dynamics, and is the most computational effort. 
An approach for thin films, described by Huang^^ is used 
here, where we need effective Green's functions that act 
in 2D on the magnetization M(x, y). This is somewhat 
different from that used in Refs. [l5| and [13, where the 
in-plane part of H^^ was calculated by first estimating 
the magnetic charge density pM- Here, it is preferred to 
calculate H^'^ directly from the field M, which has less 
steps, and is found to result in extremely precise energy 
conservation in the absence of damping. 

By applying an integration by parts, and throwing out 
a surface term outside the magnet, the solution for the 
magnetic potential is first written as an operation on M: 



$M(r) = / d'r' V'G3D(r - r') • M{r'). 



(12) 



One can notice that this involves the propagator for the 
dipole potential, that is. 



V'G3D(r-r') = 



r r 



Airlr — r' 
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(13) 



is the function whose product with a source dipole at 
position r' gives the magnetic potential at r due to that 
dipole. 

To proceed further, it is useful to consider the contribu- 
tions to the vertical (z) and horizontal (xy) components 
of separately. Consider a source cell centered at 
{x',y'), and the vertical component of it generates, 
due to M'^ = Mz{x',y'), at an observer position {x,y). 
The usual procedure is to sum over the source point z' 



and average over the observer position z. One has the 
contribution from this cell, of area dA' — dx' dy' , 



d<^>M = 



dA'M', 
47r 



dz' 



iz-z') 
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P + {z- z'Yfl"^ 



(14) 



where 8 = L/2 and the notation = (a; — a;')^ + (T/ — j/')^ 
is used. The integration gives 



M 



-dA'M'^ 
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^/r^+{z + S)^ y/P + {Z - 



(15) 

This would also be obtained exactly the same if starting 
from the magnetic surface charge density. Then, its neg- 
ative gradient with respect to z gives the contribution to 
the demagnetization field. If we also do the averaging 
over the observer position z, these two operations undo 
each other. The field averaged in the observer cell posi- 
tion is 



{dH, 



Mzl 



\ I dz-^^M = - ^d^ 
L J_s dz J 



M 



+S 



(16) 



Evaluation of the limits, and then including a sum over 
the source point r' = {x' , y'), shows that the field is deter- 
mined by convolution with an effective Green's function 
in 2D, 



Huzir) 



Gzzir) 



d^r' G,,(r-r')M,(r') 



1 



1 
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(17) 



(18) 



In these expressions, it is understood that the positions r, 
r' and the displacement between the two, f = r — r', are 
now two-dimensional. The expression for Gzz is divergent 
at zero radius. However, it is a weak divergence that 
can be regularized for the computation on the grid, by 
averaging over the cell area. For removing the divergence 
at f = 0, averaging over a circle of area equal to the cell 
area replaces the value of Gzz(O) by the longitudinal 
demagnetization factor for cylinder of length L and 
radius Tq = a/ ^/tt. So we set 

Gzz{Q) = {Gzz)o = -N, = -1 + ro - VT^TTf 

(19) 

The "o" subscript refers to averaging over the circle of 
radius Tq- See Ref. [l^ for more details. Note that Gzz 
is always negative; it correctly gives the demagnetization 
field opposite to the magnetization M which generated 

For the in-plane components of i/*^, a similar proce- 
dure can be followed. Due to symmetry considerations, 
only and My can contribute. One can start by finding 
the magnetic potential. 



d^ 



dA' , jx - x')M'^ + {y~ y')M'y 
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The integration over the source vertical coordinate z' 
gives 



M 



^ [ix-x')M', + iy^y')M;^] 



z + S z — S 



+ {z + (5)2 + (z - Sy 



(21) 



The averaging over the observer point z can be carried 
out, and gives, 



(ci$M) = j ^dzd^M{z) = 



Vf2 +i2 _ If I 



X [(x-a;')M; + (2/-y')M;]rf^'- (22) 

Finally, the in-plane gradient leads to the in-plane de- 
magnetization components. Including also the z compo- 
nents, the demagnetization field averaged in the observer 
cell is obtained from 

Ht\v) = jd'v' J2 G„^(r-r')M^(r'). (23) 

I3=x,y,z 

The elements of the Green function needed here are found 
to be 



Gxy{v) 



L 



27rf4 VVr2 + L2 
L 27^2-^X2 +r 



Vf2 



xy 



Vf2 



L2 



,(24) 
(25) 



The element Gyy is obtained from G 



and y indices, and Gyx = C?a: 



by swapping x 
One can verify that these 



matrix elements go over to those for the far-field of a 
point dipole, in the limit f — > cx). 

These transverse elements of G also are not defined at 
zero radius, because an implicit assumption in the deriva- 
tion is that the observation point is outside of the source 
cell. There needs to be an internal demagnetization effect 
within a cell even for a transverse magnetization such as 
^ Q ov My ^ 0. For long thin cells with L :^ a, 
this internal transverse demagnetization factor would be 



approximately 



As a better alternative. 



we set Gxy{0) = 0, and replace G^xiO) and Gyy{0) with 
the transverse demagnetization factor of a cylinder with 
cross-sectional radius Tq = al^^ 



G..(0) = G,,(0) = -Nx = ^ (/L^ 



(26) 



In this way, the internal demagnetization components of 
the computation cells satisfy the requirement + Ny + 
Nz = 1, while making Gxx{^) and Gyy{0) consistent with 
the regularization done for Gzz(O). 

The above results show that H^^ is found by convo- 
lution of the 2D Green's operator, as a matrix, with M. 
The calculation can be made faster by using a fast Fourier 
transform (FFT) approach,?^ which replaces the convolu- 
tion in real space with multiplication in reciprocal space. 



Of course, the simplest FFT approach requires a grid 
with a size like 2" x 2", where n is an integer. Our 2D 
system is a circle of radius R — Na {N is the size in 
integer grid units). For the FFT approach to work, so 
that the system being simulated is a single copy of the 
circle with no periodic interactions with the images, one 
can choose the smallest n such that 2" > 27V. By making 
the FFT grid at least twice as large as the circle to be 
studied, the wrap-around problem, due to the periodicity 
of Fourier transforms, is avoided in the evaluation of the 
convolution. The FFT of the Green's matrix, which is 
static, is done only once at the start of the calculation. 
During every time step of the integrations, however, the 
FFT of the magnetization field components must be car- 
ried out, for every stage at which the demagnetization 
field is required. Of course, the inverse FFTs to come 
back to H^^ are needed as well in every stage of the time 
integrator. 



III. THE DYNAMICS AND UNITS 

A. Zero temperature 

The zero-temperature undamped dynamics of the sys- 
tem is determined by a torque equation, for each cell of 
the micromagnetics system. 



dfk 
~dt 



= j^li X Bi 



(27) 



Here Bi is the local magnetic induction acting on the 
j*^ cell, 7 is the electronic gyromagnetic ratio, and the 
dipole moment of the cell is fli = La^MgiJii. The local 
magnetic induction can be defined supposing an energy 
—fli ■ Bi for each dipole, with 



SH 



J 



isn 

/i 5 fli Lo?Ms 



(28) 



The sum over j contains only sites z[i) that are near- 
est neighbors of site i. This dimensionless induction bi 
used in the simulations is converted to real units by the 
following unit of magnetic induction. 



Bn = 



J 



2A 



La? Ms a^Ms 



(29) 



For computations, the dynamics is written in terms of the 
dimensionless fields, also scaling the time appropriately: 



d-rfii 



rhi X bj 



iB^t. 



(30) 



This means that the unit of time in the simulations is — 
(7Bo)~^- For Permalloy with A = li pJ/m, Ms = 860 
kA/m, one has Aox ~ 5.3 nm. In our simulations we put 
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the transverse edge of the cells as a = 2.0 nm. Then using 
the gyromagnetic ratio, 7 = e/rrie ~ 1.76 x 10^^ s~^, 
the computation units are based on fioMg — 1.08 T and 
Bq « 7.59 T. This large value for Bq is the scale of the 
local magnetic induction due to the exchange interaction 
between the cells. The time unit is then to « 0.75 ps; 
a frequency unit is fo = 7-B0 = 1.336 THz. We may 
display frequency results, however, in units of ^jMs ~ 
15.1 GHz for Permalloy, as this expression is equivalent 
to in CGS units. For the disk sizes used here, typical 
periods of the vortex gyrotropic motion are around tq ~ 
4000, which then corresponds to dimensionless frequency 
ly = 1/tg ~ 2.5 X 10~^, and hence, physical frequency 
f = iyfo- 0.3 GHz. 

In some cases we also need to include Landau-Gilbert 
damping, with some dimensionless strength a. Then this 
is included into the dynamics with the usual modifica- 
tion. 



— — — rrii X bi — ami x ymi x bij 



(31) 



The zero temperature dynamics was integrated numer- 
ically for this equation, using a standard fourth-order 
Runge-Kutta (RK4) scheme. Typically, a time step of 
At — 0.04 was found sufficient to insure the correct en- 
ergy conserving dynamics (when a = 0) and result in 
total energy conserved to better than 12 digits of preci- 
sion over 5.0 x 10^ time steps in a system with as many as 
4000 cells. To get this high precision, however, it is nec- 
essary to always evaluate the full demagnetization field 
at all four intermediate stages of the individual Runge- 
Kutta time steps. 



B. Finite temperature: Langevin dynamics 

For non-zero temperature, the dynamics is investigated 
here using a Langevin approach. This requires including 
both a damping term and a stochastic torque in the dy- 
namics; together they represent the interaction with a 
heat bath. The size of the stochastic torques is related 
to the temperature and the damping constant, such that 
the system reaches thermal equilibrium. 

It is reasonable to think of the dynamics depending on 
stochastic magnetic inductions bs, in addition to the de- 
terministic fields bi from the Hamiltonian dynamics. For 
the discussion here, suppose we consider the dynamics 
of one computation cell, and suppress the i index. The 
dynamical equation for that cell's rh, including both the 
deterministic and random fields, is 



dm 
d7 



b + b. 



m X (6 + 6, 



(32) 



The first term is the free motion and the second term is 
the damping. Alternatively, the dynamics can be viewed 
as that due to the superposition of the deterministic ef- 
fects (due to b) and stochastic effects (due to bg). 



For a given temperature T, the stochastic fields estab- 
lish thermal equilibrium, provided the time correlations 
satisfy the fiuctuation-dissipation (FD) theorem. 



{b^{T)b^'{T'))^2arS,y 5{t-t'). 



(33) 



(5a A' is the Kronecker delta and the indices A, A' refer to 
any of the Cartesian coordinates; 6{t — r') is a Dirac 
delta function. The dimensionless temperature T is the 
thermal energy scaled by the energy unit J, 



T - 



kT 



kT 
2AL' 



(34) 



where k is Boltzmann's constant. The fiuctuation- 
dissipation theorem expresses how the power in the ther- 
mal fiuctuations is carried in the random magnetic fields. 
In terms of the physical units, the relation is 

-ffi{B^{t)B^'{t')) = 2akTSxx' 6{t - t'). (35) 

where /i = La^Mg is the magnetic dipole moment per 
computation cell. 



Time evolution with second order Heun (H2) 
method 



The Langevin equation p2l) is a first-order differential 
equation that is linear in multiplicative noise, li y — 
y{T) represents the full state of the system (a vector of 
dimension ZN , where TV is the number of cells), then the 
dynamics follows an equation of the form 



dT 



f[r.v{r)]+Is[r,y{T)]-b,{T). 



(36) 



The vector function / is the deterministic time deriva- 
tive and the vector function fa determines the stochastic 
dynamics; 6^ represents the whole stochastic field of the 
system. An efficient method for integrating this type of 
equation forward in time is the second order Heun (H2) 
method. ^'^'^^ That is in the family of predictor-corrector 
schemes and is rather stable. It involves an Euler step as 
the predictor stage, and a corrector stage that is equiv- 
alent to the trapezoid rule. Some details of the method 
are summarized here, to indicate how the stochastic fields 
are included, and to show why it is used rather than the 
fourth order Runge-Kutta method (the latter seems dif- 
ficult to adapt to the stochastic fields). 

We use the notation j/„ = 2/(t„) to show the values at 
times r„ = rtAr, according to the choice of some inte- 
gration time step Ar. Integrating Eq. (I36p over one time 
step gives the Euler predictor estimate for ?/(t„ -|- At): 

Vn+i = Vn + /(T„,?;„)Ar-|-/s(r„,y„) • {crsW^)- (37) 

The last factor, UsWm is introduced to represent the time- 
integral of the stochastic magnetic inductions. Ug is a 
variance and w„ represents a vector of 3A'^ random num- 
bers, one for each Cartesian component at each site of the 
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grid. Consider, say, the result of integrating the equation 
of motion for just one component for one site: 

jy\rhi{T)~^a,wl. (38) 

The physical variance CTs needed for this to work cor- 
rectly, must be determined by the FD theorem. For this 
individual component at one site, the squared variance is 

dr j dr' {bUrKir')). (39) 

Now applying the FD theorem to this gives the required 
variance of the random fields, that depends on the time 
step being used: 

as = %/2ar At. (40) 

This means that individual stochastic field components 
6g (r), integrated over one time step, are replaced by ran- 
dom numbers of zero mean with variance as, as used 
above. 

For the corrector stage, the points ?/„ and yn+i are used 
to get better estimates of the slope of the solution. Then 
their average is used in the trapezoid corrector stage: 

Vn+l = yn + ^ 2/„) + /(t„+i, yn + i)] At (41) 

+ ^ [/s(T„,y„) + fs{Tn+l,yn+l)] ' {cTsWn)- 

The error is of order 0{{At)^), hence it is a second or- 
der scheme. Note that the same vector of 3iV random 
numbers Wn used in the predictor stage are re-used in 
the corrector stage, because it is the evolution over the 
same time interval. 

In the coding for computations, one does not use the 
explicit form of the functions / and fs. Rather, at each 
cell, first one can calculate the deterministic effective field 
bi based on the present state of the system. Its effect in 
the dynamics will be actually proportional to its product 
with the time step, i.e., it gives a contribution Am^ cx 
biAr. Of course, the stochastic change in this same site 
will be proportional to the stochastic effective field, which 
is some a^Wi for that site, where Wi = (wf , wf , wf) . So 
the total change at this site is linearly determined by a 
combination, 

A■m^ cx cji, cji EE fejAr -f- agWi. (42) 

An effective field combination gi acts in this way both 
during the predictor and the corrector stages. In either 
stage, a dynamic change in a site is given by a simple 
relation, 

Arhi = x [g^ - a{rhi x gi)] . (43) 



Of course, the predictor stage uses the last configuration 
of the whole system to determine all the bi, while the 
corrector finds the needed bi based on the predicted po- 
sitions. And, the corrector actually does the average of 
Ami from the Euler stage and the second estimate from 
the corrector stage. The same random numbers w„ used 
in the predictor stage are used again in the corrector, for 
a chosen time step. 

The integration requires a long sequence of quasi- 
random numbers Wn ■ It is important that the simulation 
time does not surpass the period of the random num- 
bers. We used the generator mzranlS due to Marsaglia 
and Zaman,^^ implemented in the C-language for long 
integers. This generator is very simple and fast and has 
a period of about 2^^^, and is based on a combination of 
two separate generators with periods of 2^^^ and 2^^ . 



IV. VORTEX STATE PROPERTIES AND 
ZERO-TEMPERATURE DYNAMICS 

The dynamics at zero temperature, calculated with 
RK4, was used to check basic vortex dynamic proper- 
ties such as the stability and gyrotropic mode frequency. 
We also used the Langevin dynamics calculated with sec- 
ond order Heun method to include finite temperature to 
see the primary thermal effects for some specific vortex 
initial configurations. For some of these studies, it is ex- 
tremely beneficial to produce a well-formed initial vortex 
state in some desired location without the presence of 
spin waves. 

An initial vortex state is prepared first in a planar con- 
figuration of positive vorticity q — +1, namely, in-plane 
magnetization angle </> = tan^^ my /nix given by 

Mx,y) = gtan"^ x - ^ 

y-ya 

(The negative vorticity state g = — 1 is destabilized by 
the demagnetization field, so there is no reason to con- 
sider it.) This is the profile of a vortex centered at 
position {xo,yo). The out-of-plane component here is 
ruz — 0, however, the stable vortex state has a nonzero 
out-of-plane component close to rriz = ±1 at the vor- 
tex core (polarization p — ±1). This stable vortex state 
was reached by the local spin alignment procedure^- for 
a vortex at the constrained position {xQ,yo), described 
in Ref. ,15. Briefiy, that is a procedure where each rhi 
is aligned along its local induction bi, and the process 
is iterated until convergence. The constraint is applied 
as extra fictitious fields included with the Lagrange mul- 
tiplier technique, that force the desired vortex starting 
position. This procedure helps to remove any spin waves 
that would otherwise be generated starting from any ar- 
bitrary initial state. This state would be a perfect static 
state if generated in the center of the disk. When gen- 
erated off-center, the dynamics associated with its mo- 
tion still is able to produce some spin waves. A cleaner 
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vortex motion can be generated if there is a weak damp- 
ing applied {a ~ 0.02) over some initial time interval 
(r 1000). After that, the system can be let to evolve 
in energy-conserving dynamics, if needed. 

This relaxed vortex state develops either positive or 
negative out-of-plane component, including some small 
randomness in the initial state before the relaxation. If 
rriz ~ +1 (—1) in the vortex core region, the vortex has 
positive (negative) polarization and a positive (negative) 
gyrovector G = G^, defined from 



G = 2t:Q — 2, Q = qp. 

7 



(45) 



7 is the electron gyromagnetic ratio and toq — A'/a^ = 
LMs is the magnetic dipole moment per unit area. The 
integer Q = ±1 defines the quantized topological charge 
that determines the two allowed discrete values of the gy- 
rovector. To a good degree of precision, the vortex states 
studied here obey a dynamics for the vortex velocity V 
described by a Thiele equatiouj^iii ignoring any intrinsic 
vortex massii or damping effects. 



F + Gx V = 0. 



(46) 



This equation comes from an analysis of the Hamilto- 
nian dynamics of a magnetic system^^i^ in which the 
vortex excitation profile preserves its shape but moves 
with some collective coordinate center position X(i), 
with Y{t) = X(i). The force F is the gradient of the 
potential experienced by the vortex. The force points 
towards the nanodisk center, and can be approximated 
by some harmonic potential with force constant kp, for 
a vortex at distance r from the center, 



F = -k 



pr r. 



(47) 



Hence, the presence of the gyrovector leads to the well- 
known gyrotropic (or uniform circular) motion. Solving 
for the vortex velocity results in 



V 



G 



^kpT 
2ttQLMs 



(48) 



G includes the sign of the gyrovector (vector G points 
perpendicular to the plane of the disk, and it has only a 
z component). Thus, the vortices generated with positive 
(negative) gyrovector move clockwise (counterclockwise) 
in the xy plane. Furthermore, the angular frequency of 
this gyrotropic motion is given by a related equation. 



LOG 



V 



■ykf 



kp 

~G ^ 2tiQLMs 



(49) 



The force constant has been estimated theoretically from 
the rigid vortex approximation^ and from the two- vortex 
modelii Below, we determine kp numerically from re- 
laxed vortex states^ (a flexible vortex). The frequency 
in Eq. (|49|) applies to the stable vortex states. If the 
disk is too thin, the vortex could be unstable; this pro- 
duces an outward force F, and results in the gyrotropic 
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FIG. 1: Vortex motion with damping, at zero temperature. 
This is clockwise motion for a vortex with positive {-\-z) gy- 
rovector, starting from the dot on the a;-axis. The vortex 
performs gyrotropic motion of decreasing radius and increas- 
ing frequency as it moves towards the disk center, r = (0, 0). 



motion in the "wrong" direction. Thus it is easy to iden- 
tify whether a vortex is stable or unstable from a short 
integration of its dynamics. 

In the time and frequency units applied in the sim- 
ulations, the dimensionless gyrotropic frequency VLg is 
obtained from 



— ^cto — 



LOG kpa^ 
^ ^ AttLAQ ' 



(50) 



The negative sign shows that vortices with a negative 
gyrovector {Q = —1) have a counterclockwise rotational 
motion; the opposite sense holds for positive gyrovector. 
The force constant kp increases with thickness L but de- 
creases with disk radius R. Therefore, in the simulation 
time units, the gyrotropic frequency could depend pri- 
marily on their ratio, L/R. 

For detection of the vortex motion, one method is to 
measure the spatially averaged magnetization. 



N ^ 



(51) 



This is a useful measure of vortex gyrotropic motion, es- 
pecially for experiments, where it may not be possible to 
observe the rapidly changing instantaneous vortex core 
position. However, (m) can show rotational oscillations 
even when no vortex is present. Thus, we need instead a 
measure of the vortex core position based on the location 
of the vorticity charge center. 

The vorticity center position r„ is the point around 
which the in-plane magnetization components give a di- 
vergent curl. That is, a continuum magnetization field of 
a vortex located at position r„, with in-plane angle 0(r), 
would be expected to have the curl, 



V X V0(r) = 2-Kzi{v - r„). 



(52) 
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FIG. 2: (Color online) For the vortex motion in Figure [T] 
the phase relationship between perpendicular components of 
position and in-plane magnetization. 



When used on the discrete grid of cells, the vorticity cen- 
ter falls between the four nearest neighbor grid cells that 
have a net 2tt circulation in cf). However, this discretely 
defined position always jumps in increments of the cell 
size a, hence, it cannot be used directly. Instead, we use 
an average position weighted by the squared mf compo- 
nents, of only those cells near the vorticity center: 



E 



\<iX, 



,(mf) 



z \2 , 



(53) 



The Ti are the cell positions and the sum is restricted to 
those cells within four exchange lengths of the vorticity 
center. The center of the nanodisk is the origin, (x, y) — 
(0,0). Including this cutoff in the sums helps to reduce 
the contributions from other oscillations in the system 
(i.e., spin waves) that are not directly associated with the 
vortex position. By weighting with (to|)^, the position Tc 
is able to change smoothly as the vortex moves, especially 
at T = 0, in contrast to the discrete vorticity center r^. 
It is a reasonable estimate of the mean location of out- 
of-plane magnetization energy of the vortex, i.e., close to 
the vortex core position. The rriz-weighted position Vc 
and the vorticity center r^, are usually within one lattice 
constant. This measure is supplemented by observing the 
actual magnetization field when there is any doubt about 
the presence or stability of the vortex. 



A. Gyrotropic frequencies in circular disks 

Calculations were carried out for circular disks of thick- 
ness 5.0 nm, 10 nm and 20 nm (L = 2.5a, 5a, 10a, all with 
a = 2.0 nm) for radii 30 nm, 60 nm, 90 nm and 120 nm. 
The stability of the vortex state is easily checked for a 
given geometry, by starting from a relaxed vortex at some 
radius near half the radius of the disk. Including a weak 



FIG. 3: (Color online) Typical motions of the vortex core 
coordinate Xc{t) at zero temperature, for circular disks of 
thickness L = 10 nm with different radii (shifted vertically 
from Xc = for clarity). The damping a = 0.02 was turned 
off at time r = 1000. Periods were calculated from the energy- 
conserving motion after t > 1000. The motion of t/c(T) is 
similar but shifted a quarter of a period. 



damping a = 0.02, it is necessary only to run a short sim- 
ulation of the dynamics and observe whether the vortex 
moves in the direction given by the Thiele equation;^ 
Eq. (gSD. 

For example, with i? = 30 nm, L — 5.0 nm, a vor- 
tex was initially relaxed at a position {xq^i/q) = (16,0) 
nm, and then the dynamics was started, including damp- 
ing a ~ 0.02 in the RK4 method. In this case the vor- 
tex is very stable and spirals into the center of the disk, 
see Figures [1] and [2l The instantaneous vortex displace- 
ment on one axis, scaled by disk radius, takes approxi- 
mately the same magnitude as the perpendicular in-plane 
component of {m), such as Xc/R and (my) in Figure [21 
Another feature is that the period of rotation becomes 
less as the vortex moves inward. The first few peri- 
ods are At = 6000, 3140, 2580, but the later revolutions 
have an average period tq ~ 2020 (1.51 ns, frequency 
/g = I/tg = 0.661 GHz for Py). 

Other similar dynamics calculations were done at vari- 
ous disk sizes, but turning off the damping a — 0.02 after 
T — 1000, see Figure [31 This initial damped motion is 
used to remove spin waves that might be generated when 
the vortex is initially released, after being relaxed at a 
desired starting position. Once the damping is turned 
off, the dynamics is energy conserving. Because we are 
later interested in small movements near the disk center, 
the initial position was taken as {xq, yo) = (2a, 0), using a 
lattice constant a = 2.0 nm. These simulations result in 
very smooth circular motion of the vortex center v,. (Fig. 
[3]), from which very precise estimates of the gyrotropic 
period tq were determined by following the motion for 
typically five to ten periods. The resulting frequencies 
/g, in units of ^7Ms, are shown versus aspect ratio 
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B. Relation to force constant kp 
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FIG. 4: (Color online) Zero-temperature vortex gyrotropic 
frequency fc for various disk radii R, versus aspect ratio L/R. 
[For Permalloy, ^^yMs ~ 15.1 GHz]. The computation cell 
size is a = 2.0 nm. The vortex state is unstable below a 
minimum disk thickness, as expected due to the diminished 
restoring forces from the reduced edge area. The dashed line 
shows the result [Eq. (I59p ] from using the linear approxima- 
tion in Eq. (f55)) for kp- 
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FIG. 5: (Color online) Vortex force constant kp scaled by 
disk thickness, versus disk aspect ratio. These were obtained 
by assuming a parabolic potential for vortex motion within 
the disk. The dashed line indicates that the slope of this 
relationship is close to 1/4 for some range of parameters, Eq. 
(|55[) . for disks of adequate thickness. Cell edge is a = 2.0 nm. 



L /R in Figure S) The scale is also given there for the 
parameters of Permalloy, for which ^jMg « 15.1 GHz. 
One can note the obvious feature, that the gyrotropic fre- 
quency goes to zero at some minimum thickness needed 
for vortex stability. 



The vortex restoring force constants kp were estimated 
based only on static energy considerations. We compared 
the total system energy with the displaced vortex, U{x), 
taking x — 2a, with the energy for the vortex at the disk 
center, U{0). It is known that the vortex potential is 
close to parabolic, as long as the vortex displacement is 
small compared to the disk radius The force constant 
is then estimated simply by solving 



U{x) = U{0) + 



kp x^. 



(54) 



The energies applied in this equation are those obtained 
after the vortex is relaxed by the Lagrange-constrained 
method. These calculations are relatively fast because 
there is no need to run the dynamics. The raw force 
constants were obtained for a wide variety of disk sizes. 
Generally, we find that kp increases faster than linearly 
with disk thickness L and decreases with disk radius R. 

It is expected that the force constant should scale 
somewhat with the aspect ratio, LjU. Further, the 
Thiele equation suggests that the ratio kp/L is most rel- 
evant in determining ujq [see Eq. (|49p ]. Therefore, we 
show kp/L versus L/R in Figure O which presents a re- 
lationship somewhat close to linear, with a slope near 
1/4. Thus we can write as a rough approximation (far 
enough from the critical disk thickness for vortex stabil- 

ity), 



The last form, obtained by applying the definition of ex- 
change length, is preferred because the vortex restoring 
force ultimately is due to the demagnetization fields gen- 
erated by Ms. 

One can check whether these force constants are consis- 
tent with the gyrotropic frequencies found in the dynam- 
ics. If the Thiele equation applies to this motion, then 
the gyrotropic frequencies must be linearly proportional 
to kp/L, [Equations and ([50)) ]. Therefore we have 
plotted the dimensionless frequency Q.g versus kp/L in 
Figure [S) For the wide variety of disk sizes studied, all 
points in this plot fall on a single line of unit slope, ex- 
actly consistent with the Thiele equation. This shows 
that the calculations of the dynamics over fairly long 
times (many periods) are completely consistent with the 
force constants found only from static energy consider- 
ations. It further implies that we can safely use static 
energy calculations to predict dynamic properties. This 
is based on the assumption of an isotropic parabolic po- 
tential in which the vortex moves. There may be some 
limitation to this idea, however, only because the poten- 
tial will deviate from parabolic for larger displacements 
from the disk center. 

These results are consistent with the two-vortices 
model applied by Guslienko et al."^ With the boundary 
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FIG. 6: (Color online) The dimensionless gyrotropic frequen- 
cies (found from dynamics) versus force constant scaled by 
disk thickness. The dashed line of unit slope is Eq. (|50p . This 
verifies the dynamics of the Thiele equation, and shows the 
complete consistency between the static energetics and the 
dynamics. Cell edge is a = 2.0 nm. 



parameter ^ = 2/3 and the initial susceptibility at small 
aspect ratio being x(0)~"'^ w 9.98L/R, their result (con- 
verted to SI units by factor f^) is approximately 



At: 



(56) 



Our results have a somewhat weaker potential, which is 
to be expected because the numerical simulations allow 
for a wider range of possible deformations of the vortex 
structure than is possible in an analytic approximation. 
In addition, our numerical results include the destabiliza- 
tion of the vortex at sufficiently small L/R, hence, it is 
impossible to fit any straight line for kp/L vs. L/R down 
to arbitrarily small aspect ratio, see Figure [5l 

We showed above that the gyrotropic frequencies z/g 
are exactly linearly proportional to kp/L, hence, this im- 
plies that the frequencies also scale close to linearly with 
L/R. Combining our fit of kp with relation ([50]) then 
shows that roughly, the dimensionless angular frequency 
magnitude is 



1 L L 
Idtt R R 



In physical units, this is 



= -/Bona w O.UOjfioMs^. 



Then the frequency comes out 
fa 



L 



^« 0.280 (^jMA- 
27r V47r ' / R 



(57) 



(58) 



(59) 



The dashed line in Figure U shows Eq. ((5^ compared 
with data from various disk sizes. These frequencies are 



smaller than those in the rigid vortex model;^ and only 
slightly smaller than those for the two-vortices modelfi 
However, this result fits quite well with the experimental 
data presented in Ref. ^ by also using the higher value 
for the gyromagnctic ratio, 7 = 1.85 x 10^^ s~^ T"-'^, 
in conjunction with saturation magnetization still at the 
value Ms = 860 kA/m. The calculation here can be 
considered as that for a more fiexible vortex. The mag- 
netization at the edge of the disk adjusts itself to try to 
follow the boundary. The magnetization can also adjust 
itself, to a lesser extent, in the vortex core region. These 
effects lead to lower force constants and therefore lower 
gyrotropic frequencies. 

These results show that the adapted 2D methods ap- 
plied here give reliable results, consistent with experi- 
ment and with the two- vortices analytic calculation of the 
gyrotropic frequencies. We note that the smaller value of 
cell constant used here (a = 2.0 nm) is important for the 
simulation to correctly describe the magnetization dy- 
namics in the vortex core. Of course, this then imposes 
a limitation on the system size that can be studied. 

These results confirm the basic dynamic properties, 
that the vortex resonance frequency ujg diminishes with 
increasing dot radius, and increases with increasing dot 
thickness. A wider dot has a weaker spring constant kp 
in its potential, U{r) = U{0) + ^kpr^ , leading to the re- 
duction of its resonance frequency. Similarly, in a thicker 
dot, the greater area at the edge produces a larger restor- 
ing force, leading to a higher resonance frequency. 



V. THERMAL EFFECTS IN VORTEX 
DYNAMICS IN CIRCULAR DISKS 

In the following part, the effects of thermal fluctuations 
on the vortex dynamics are considered. We consider two 
basic situations left to evolve in time via Langevin dy- 
namics: (1) A vortex started off-center, and (2) a vortex 
started at the minimum energy position, the center of 
the disk. In the latter case, the question is whether ther- 
mal fluctuations alone are sufficient to initiate gyrotropic 
motion. If so, we can also study its frequency and range 
of motion. In all simulations we used cell size a = 2.0 nm 
and damping parameter a = 0.02 . 



A. Vortex initially off-center 

For the same system used above [R = 30 nm, L = 5.0 
nm], the same initial condition was used, with vortex 
at (xqjJ/o) = (16,0) nm, but a finite temperature cor- 
responding to Permalloy at 300 K was considered. The 
dynamics was solved now by the H2 scheme. The scaled 
temperature depends on the thickness L of the disk and 
the exchange stiffness A of the material. The energy 
unit here is J = 2AL = 130 zJ, while 300 K corre- 
sponds to kT = 4.14 zJ, so the scaled temperature is 
T = kT/ J — 0.032. The x-component of the vortex po- 
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FIG. 7: Vortex motion in Py at room temperature (300 K), 
starting from an initial displacement of 16 nm from the disk 
center. The y-component of average magnetization in the disk 
is correlated to the a;-component of the vortex position. 



FIG. 8: Spontaneous gyrotropic vortex motion in Py due to 
thermal fluctuations at 300 K, starting from a vortex at the 
center of the disk. 
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FIG. 9: Thermal power spectrum of the in-plane magnetiza- 
tion fluctuations due to spontaneous gyrotropic vortex motion 
in Py at 300 K, for the motion in Figure [S] 



sition versus time is shown in Figure [71 In this case, the 
vortex still spirals tovifards the center of the disk, how- 
ever, thermal fluctuations remain present in the motion 
even at time r = 60000 (« 90 ns), 25 revolutions later. 
The range of the motion there remains close to ±6 nm. 
The time dependence of (my) (most closely related to Xc) 
is also shown in Figure [3 it also shows an effect persisting 
at the 25% level out to t = 60000. Note that at zero tem- 
perature, the time-scale for relaxation (Figure [2]) was on 
the order of r ~ 20000. This shows that thermal forces 
apparently are able to maintain the gyrotropic motion 
to very long times. The average period of the motion is 
TG ~ 2278 (1.705 ns, frequency / = 1/tg ^ 0.586 GHz 
for Py), showing that the temperature also softened the 
potential experienced by the vortex. 



B. Vortex initially at disk center 

The same system is used [R = 30 nm, L = 5.0 nm], 
but this time the vortex was initiated at the center of 
the disk, {xo,yo) = (0,0). At zero temperature, such an 
initial state is static. Instead, the dynamics correspond- 
ing to Py at 300 K was considered (scaled temperature 
T = 0.032). Any thermal fluctuations can move the vor- 
tex core off-center, and if that happens, gyrotropic mo- 
tion can initiate spontaneously. This indeed happens, 
as can be seen in the vortex core position rc(T) plotted 
in Figure [S] It needs to be stressed that these vortex 
motions of the order of ±4 nm, and magnetization fluc- 
tuations on the order of ±15%, occur without the ap- 
plication of any external magnetic field. The motion is 
sufficiently coherent that it can be followed for dozens 
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R=120 nm, L=20 nm, 300K (kT/J=0.008) 




FIG. 10: Spontaneous gyrotropic vortex motion, due to ther- 
mal fluctuations, in a 20 nm thick Py disk at 300 K, with 
the vortex starting at the center of the disk. The natural pe- 
riodic motion executes 32 revolutions in this time sequence, 
with period tg ~ 1870. 



of rotations. The gyrotropic motion was followed out to 
twice the time shown in the plots. An average over 24 
rotations results in a period tq = 2250, corresponding to 
1.68 ns or a frequency / = 0.594 GHz. To verify this, we 
also show the power spectrum of the in-plane magnetiza- 
tion oscillations in Figure El This was obtained by taking 
time FFTs of {mxir)) of length 256 points at different 
starting times in the data out to r = 120000 and aver- 
aging their absolute squares. The middle peak in Figure 
ini falls at dimensionless frequency ly w 4.52 x lO"**, cor- 
responding to physical frequency / — v/to = 0.600 GHz, 
consistent with the estimate from counting oscillations. 
There is some structure in the FFT, possibly the beating 
between three different primary frequencies, that causes 
the amplitude of the oscillations to wax and wane. 

The spontaneous gyrotropic vortex motion takes place 
for a wide range of system sizes that were tested. Another 
example is given for a larger system [R = 120 nm, L = 20 
nm] in Figure I10[ where the vortex core displacement is 
displayed. An interesting feature is apparent. The gy- 
rotropic motion loses its phase coherence at times, lead- 
ing randomly to brief intervals of dramatically changed 
amplitude. This is only one example; in other time se- 
quences for other system sizes, this behavior is particu- 
larly intermittent and random. For the same simulation. 
Figure [TT] also shows both components of vortex core po- 
sition and both components of the average in-plane mag- 
netization, zoomed in to show details at earlier times. 
Here one can see the quarter-period phase difference be- 
tween X and y components for the vortex position as well 
as for the magnetization. In addition, the magnetiza- 
tion exhibits a high-frequency oscillation with a period 
of about At « 125 on top of the gyrotropic oscillations. 
This can be expected to be spin wave excitations that 
are excited thermally together with the vortex gyrotropic 
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FIG. 11: (Color online) For the spontaneous gyrotropic vortex 
motion in Figure [TOl [R = 120 nm, L = 20 nm Py disk at 300 
K], details of the motion at earlier times. The vortex started 
at the center of the disk. There is a high-frequency spin wave 
oscillation apparent in the magnetization dynamics, excited 
together with the gyrotropic motion. 



motion. 

To confirm the identity of these spin wave oscillations, 
we also show in Figure[T2]the power spectrum in net mag- 
netization component rrix, from a longer simulation out 
to time T — 2.5 X 10^ . The vertical scale has been zoomed 
in to bring out the appearance of a doublet with frequen- 
cies of 9.3 GHz and 11.4 GHz, for Permalloy parameters, 
while the gyrotropic frequency is only 0.71 GHz. A spin 
wave doublet with azimuthal quantum numbers m = ±1 
(wavefunction varying as "0 gjm^ a,round the disk cen- 
ter) has been discussed in Ref. [l^. The doublet is pre- 
dicted to have a splitting'^^ of A/ = f2 — fi — 3.5/g 

and an averaged frequencj*?'^ of / = 1.8 (^^'yMg) \J^- 
For the situation here, these formulas predict A/ = 2.5 
GHz and / = 11.1 GHz, while the observed doublet has 
A/ = 2.1 GHz and / = 10.3 GHz. Although slightly 
softer, these are of the right orders of magnitude and are 
consistent with the the theoretical prediction for this dou- 
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FIG. 12: (Color online) The thermally averaged power spec- 
trum in one component of the magnetization (squared FFT) 
for the vortex motion in Figure 1101 The low frequency gy- 
rotropic mode dominates strongly over a much weaker dou- 
blet at high frequency. For Permalloy parameters (/ — 
1336GHz X v), the gyrotropic frequency is fa = 0.71 GHz 
while the components of the doublet lie at /i — 9.3 GHz and 
/2 = 11.4 GHz. 



blet. This lowest doublet relates to the presence of spin 
waves propagating azimuthally around the disk, in the 
presence of the vortex. The splitting can be attributed 
to the breaking of symmetry for the two directions of 
propagation, due to the presence of the out-of-plane mag- 
netization at the vortex core. Based on these results and 
results at other disk sizes, we then note that the primary 
deviation from a smooth gyrotropic motion is due to the 
thermal excitation of this doublet on top of the vortex 
magnetization. 



FIG. 13: (Color online) Typical spontaneous fluctuations of 
the vortex core a;-coordinate for 30 nm radius Py disks with 
various thicknesses, at 300 K. The vortex was initiated at the 
disk center. Curves are shifted vertically from Xc = for 
clarity. 
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C. Analysis of thermal vortex motion in circular 
nanodisks 

The spontaneous vortex motion at 300 K takes place 
without the application of any externally generated mag- 
netic field. Only the thermal energy is responsible for the 
motion. Indeed, both the frequency and amplitude of this 
spontaneous gyrotropic motion is determined directly by 
the temperature. Here we give some analysis and sug- 
gest where this motion might be most easily observed 
experimentally. 

For some smaller disks with i? = 30 nm, and for some 
larger disks, with R = 120 nm, Figures [13] and [14] ex- 
hibit the typical time dependence of the vortex coordi- 
nate Xc{t), for Permalloy systems at 300 K. The vortex 
was initially relaxed at the center of the disk (x = y = 0). 
As seen for the systems studied above, the gyrotropic 
motion is spontaneous, and furthermore, takes place at 
a lower frequency for thinner disks. In addition, there 
is a dependence of the amplitude of the motion on the 



FIG. 14: (Color online) Typical spontaneous fluctuations of 
the vortex core ^-coordinate for 120 nm radius Py disks with 
various thicknesses, at 300 K. The vortex was initiated at the 
disk center. Curves are shifted vertically from Xc = for 
clarity. 



disk thickness. The amplitude is observed to be larger 
for thinner disks. Also it is apparent that generally the 
amplitude is larger for the larger radius disks. This is 
somewhat difficult to analyze precisely, due to the lim- 
ited time sequences that can be obtained during a rea- 
sonable computation time. However, from knowledge of 
the force constants kp and their dependence on the disk 
geometry, the RMS range of the vortex core motion can 
be predicted. 

The statistical mechanics of the vortex core position 
X = {X{t),Y{t)) and velocity V = X can be obtained 
from the effective Hamiltonian associated with the Thiele 
equation. The Thiele equation is mathematically equiv- 
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alent to the equation of motion for a massless charge 
e in a uniform magnetic field B, with eB = — G, and 
also affected by some other force F. We can start from 
a Lagrangian that leads to the Thiele equation, using 
the symmetric gauge for the effective vector potential, 
and including a circularly symmetric parabolic potential 
(harmonic approximation) , 

L(X, X) = -iG(Xr - YX) - \kF{X^ + Y^). (60) 

The first term on the RHS is equivalent to eV • A, with 
vector potential A = ^BxXina magnetic problem; 
there is no usual kinetic energy term like ^mV^, be- 
cause the intrinsic mass is considered zero here. Only 
the z-component of the gyrovector is present, G = Gz — 
2TrpqmQj~^. Then the components of the Thiele equa- 
tion are recovered from the Euler-Lagrange variations. 



dL d dL 

dX ~ dt'dX 
dL d dL 
W ^ didY 



-kpX - GF = 0, 
-kpY + GX = Q. 



The Lagrangian is written equivalently as 



L(X,V) = -i(G X X) • V 



This leads to the canonical momentum, 



§ = -iGxX = (fy,-fx). 



(61) 
(62) 

(63) 

(64) 



This allows the transformation to the collective coordi- 
nate Hamiltonian, iJ(X, P). Following the usual pre- 
scription, we have 



1. 



1. 



H(X,P) = P-X-L = -kp^^ = -kp {X 



(65) 



Note that the derivation of the Hamiltonian does not 
depend on the choice of the gauge for the gyrovector (i.e., 
for its effective magnetic field). In Ref.'^'*, it is shown that 
the Landau gauge leads to the same result for but 
where P — GY is found to be the momentum conjugate 
to X. 

Technically this is all that is needed to analyze the 
statistics of the vortex position. By being purely po- 
tential energy, however, this Hamiltonian needs careful 
treatment. Its variation via the Hamiltonian equations 
of motion does not lead back to the correct dynamics, 
i.e., it does not give the Thiele equation. One can see 
that the difficulty is due to the fact that the position and 
canonical momentum coordinates are redundant, since 



= \GY and Py = 



-\GX. Even so, all of these 



should be considered linearly independent mechanical co- 
ordinates, and all should appear in H to give the correct 
dynamics (gyrotropic motion does not conserve X nor P, 
so both should appear in H). For that to work out, iJ 
must be expressed so that there are both potential and 
kinetic energy terms. (A similar care is needed even in 



the Landau gauge, where GY must be identified by and 
replaced as the momentum P conjugate to X?) We can 
split out half of the potential energy and redefine it in 
terms of P^ as a kinetic energy. 



7J(X,P) = lfc;.X2 



2P 

~G 



(66) 



One can easily demonstrate that the correct dynamic 
equations result only by allocating exactly half of the en- 
ergy as kinetic energy and half as potential energy. This 
then leads to the Hamilton dynamic equations for oscil- 
lations along the two perpendicular axes. For example, 
along X there is 



X = 



P. = - 



dH _ 2kFPx 
OH 



dX 



G2 ' 
= -hkpX. 



(67) 
(68) 



These give a second order equation for simple harmonic 
motion (SHO), 



X 



fr2 



(69) 



The other variations with respect to Y and Py lead to 
the same dynamics for Y. However, note that the Thiele 
equation is recovered from these dynamics only by in- 
cluding the connection (|64p that defines the canonical 
momentum in terms of the position. 

It is clear that the Hamiltonian ((66|) is the same as that 
for a two-dimensional simple harmonic oscillator with co- 
ordinate X and momentum P. For that oscillator, the 
effective spring constant is fcsHO = ^kp, and the corre- 

sponding effective mass is mgHO = ^fj- It is interesting 
to see that these lead back to the natural frequency of 
gyrotropic motion [or see Eq. (p^ ]. 



= WSHO 



fcsHO 



m-SHO 



kp 
~G' 



(70) 



Of course, as G is proportional to the disk thickness via 
the factor mo = LMs, and kp depends on both R and 
L, then this contains the various geometrical effects, es- 
pecially those associated with the vortex force constant. 

In consideration of the classical statistical mechanics, 
the important fact here is that the Hamiltonian (|65|) has 
a dynamics due to only two coordinates {X, Y) appear- 
ing quadratically. Although the dynamic equations for X 
and Y must come from the Hamiltonian (j66p of the equiv- 
alent 2D SHO, the phase space of the Thiele dynamics is 
more restricted, due to relation (IM)) between P and X. 
This forces the Thiele phase space to be only two dimen- 
sional; this does not depend on the choice of the gauge. 
As an example of that reduction of the phase space, ellip- 
tical motions are present for the 2D SHO, while the zero- 
temperature Thiele dynamics has only circular orbits. As 
we are considering thermal equilibrium, each independent 
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FIG. 15: (Color online) Average squared displacement of 
the vortex core from the disk center, versus reciprocal force 
constant. The points come from simulations out to time 
r — 2.5 X 10^; the solid lines are the predictions from the 
equipartition theorem, Eq. 1721 using the parameters for Py. 



quadratic coordinate receives an average thermal energy 
of ^kT. This gives the connection needed to predict the 
average RMS vortex displacement from the disk center. 
Specifically, for each vortex core coordinate. 



= hkT. 



(71) 



Then the average squared displacement of the vortex 
from the disk center should be 



2kT 
kp 



(72) 



These show that the average thermal energy in the vortex 
motion must be 



(i?(x,p)) = fcr. 



(73) 



Therefore, we can check that these relations actually hold 
in the simulations. The average squared displacement 
should be proportional to the reciprocal of the force con- 
stant, with the same proportionality factor (twice the 
temperature) when disks of different geometries are con- 
sidered. Some results for the average squared displace- 
ments versus reciprocal force constant in different geome- 
tries are given in FigureHSl The results depend on the be- 
havior of the force constant with disk geometry, showing 
the importance of static calculations for understanding 
the statistical dynamics behavior. The simulation data 
have a general trend consistent with Eq. [721 but there are 
large fluctuations due to the finite time sequences used, 
which is more of a problem for the systems with small 
kp- 

We can further substantiate the statistical behavior of 
the vortex core, by calculating the probability distribu- 
tion p{r) of its distance r = V X'^ + from the disk 




FIG. 16; (Color online) Probability distributions in Py disks 
of radius 30 nm at temperatures 300 K and 150 K, for the 
radial position r of the vortex, measured from the disk center, 
in units of the cell size, a = 2 nm. Solid curves are the 
theoretical expression (|74|l based on a Boltzmann distribution 
using the static force constants; points are from simulations 
out to time r = 2.5 x 10^ . 



center. Assuming that its position is governed by Boltz- 
mann statistics for Hamiltonian (j65p , the normalized dis- 
tribution from p{r)dr cx 27rr dr e^^^ is predicted to be 



p(r) = /3fci.re-^^'=^''' 



(74) 



where j3 — IkT)^^ is the inverse temperature. This dis- 
tribution also has some particular distinctive points that 
are relatively easy to check. For instance, the distribu- 
tion has a peak at the point of maximum probability, at 
the radius 



kT _ 
kp ~ V2 

In addition, the value of the function at this point is 

-1/2 



(75) 



Pn 



- P(''max) 



(76) 
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for different disk thicknesses. For these smaher systems, 
the agreement is quite good between the simulations and 
the theoretical expression, Eq. (fM]) . 

The distributions were also found in simulations for 
larger radius, see Figure [T7] for the distribution at i? = 
120 nm. In this case, the errors are considerably greater. 
This is due primarily to the larger gyrotropic period. 
Over the sampling time interval to r = 2.5 x 10^, there 
are less periods being sampled. The system has a some- 
what erratic behavior, in that the orbital radius of the 
vortex motion seems to switch suddenly between differ- 
ent values, as already mentioned. As a result, at this 
system size a greater time interval is needed to obtain a 
sample that could be considered in thermal equilibrium, 
with well defined averages. 

For thinner disks, the number of revolutions in the 
given time interval is lesser, which means the thinner 
disks may also require longer time sequences to give the 
same relative errors. Of course, the thinner (thicker) 
disks have a weaker (stronger) force constant, leading 
to the greater (lesser) amplitude spontaneous motions. 
This is clearly exhibited in the probability distributions. 
Although these aspects may be difScult to verify exper- 
imentally, the results do indeed point to much stronger 
spontaneous gyrotropic fluctuations for very thin mag- 
netic disks. In the cases where these motions were of 
greater amplitude, there may start to appear deviations 
from the distribution in (fM|. simply because the larger 
amplitude vortex motions cause the vortex to move out 
of the region where the potential is parabolic. 



FIG. 17: (Color online) Probability distributions for vortex 
radial position in Py disks of radius 120 nm, as explained in 
Figure [TBI 



We have found that the vortex core position satisfies 
this distribution reasonably well, while the vortex is 
undergoing the spontaneously generated gyrotropic mo- 
tion. There is a certain difficulty to verify this, be- 
cause very long time sequences (we used final time r — 
250000) are needed so that many gyrotropic revolutions 
are performed. During the motion, at times there are 
rather large fluctuations in the amplitude of the motion. 
The motion varies between time intervals of smooth gy- 
rotropic motion of large amplitude and other time inter- 
vals where the motion seems to be impeded, and is of 
much smaller amplitude. Even so, we were able to take 
these long sequences and produce histograms of the vor- 
tex radial position to compare with the predicted proba- 
bility distribution. An example for i? = 30 nm is given in 
Figure [TBI The temperatures are defined here by apply- 
ing the material parameters for Permalloy (that is, 300 
K corresponds to kT = 0.1592Aa, where the exchange 
stiffness for Py is A = 13 pJ/m and cell size a = 2.0 nm 
was used in all simulations). The data (points) are com- 
pared with the prediction of equation (|74l) (solid curves). 



VI. DISCUSSION AND CONCLUSIONS 

The calculations here give a precise description of 
the magnetostatics and dynamics for thin-film nanomag- 
nets, especially in the situations where a single vortex is 
present. The continuum problem for some finite thick- 
ness L has been mapped onto an equivalent 2D prob- 
lem, i.e., the modified micromagnetics adapted here. For 
high aspect ratios, L ^ 2R, the shape anisotropy is very 
strong, and this 2D system is a very good approximation 
of the full 3D problem, because it leads to the physical 
situation where the magnetization has little dependence 
on z and is predominately planar, except in the vortex 
core. 

At zero temperature, we have been able to test this 
approach and compare with the predictions for vortex 
gyrotropic motion based on the Thiele equation. This 
comparison is made possible here because the vortex force 
constants kp can be calculated from the energetics of a 
vortex with a constrained position. The application of 
the Lagrange undetermined multipliers techniquei^ for 
enforcing a desired static vortex position X has been es- 
sential in the determination oi kp. In addition, that re- 
laxation procedure also is of great utility for initiating a 
vortex at some radius while removing most of the initial 
spin wave like oscillations that would otherwise be gen- 
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erated when the time dynamics is started. As a result, 
we have been able to determine the zero temperature gy- 
rotropic frequencies for the motion of the vortex core, 
X(i), to fairly high precision. The confirmation of the 
applicability of the Thiele equation to the T = dynam- 
ics of vortex velocity V is impressive, as demonstrated 
in the straight line fit for gyrotropic frequency log versus 
scaled force constant kp/L in Figured This shows the 
complete consistency between the statics calculations of 
the force constants and the dynamics calculations of the 
frequencies, when interpreted via the Thiele equation. 

At larger disk radii, the gyrotropic frequency is found 
to be close to linear in the aspect ratio, L/R, see Eq.[58l 
The frequencies are also close to those found in the two- 
vortices model and micromagnetics calculations carried 
out in Ref. 7. The differences from those results may be 
due to the fact that we have used the cell parameter a half 
of what was used in Ref. 0- This is important, because 
the cell parameter should be sufficiently less than the 
exchange length for results to be reliable. Otherwise, if 
a is too large, the details of the energetics and dynamics 
in the vortex core cannot be correctly represented. 

At T > 0, the Langevin dynamics shows some surpris- 
ing behavior that was reported earlier in Ref. [l^, even 
when the vortex is initiated at the center of a nanodisk. 
The thermal fluctuations are indeed sufficiently strong to 
produce a spontaneous motion of the vortex core, with- 
out the application of any external held, which is not a 
simple random walk. Instead, the gyrotropic nature of 
the motion is still present, and in fact, persistent vor- 
tex rotation is the dominant feature of the motion. The 
thermal fluctuations can be viewed as a perturbation on 
top of the gyrotropic motion, however, it is the tempera- 
ture that determines the expected squared radius of the 
orbit. The orbital radius is very well described from the 
statistical mechanics of the vortex collective coordinate 
Hamiltonian (1551) . that possesses only the potential en- 
ergy associated with the vortex force constant. 

Integrations of the dynamics over very long times 
(equivalent to hundreds of vortex revolutions) shows that 
the statistics of the vortex position follows the simple 
Boltzmann distribution in Eq. (|74p . The average squared 
vortex displacement from the origin, r^^^, scales linearly 
in the temperature divided only by the force constant 
kp- This is in contrast to the vortex gyrotropic frequen- 
cies, which depend on kp/L. Thus, the results for force 
constant indirectly predict the expected position fluctu- 
ations. However, very long time sequences are needed to 
see this average behavior; over some short time intervals 
there can be large variations in the instantaneous vortex 
orbital radius. The largest spontaneous vortex position 
fluctuations will be possible in thin dots of larger radius, 
where the force constants are weakest. Even so, this is a 
small effect (RMS radii on the order of several nanome- 



ters), and it may be difficult to observe experimentally. 
As an example based only on the calculated force con- 
stants, a magnetic dot of radius R = 180 nm and thick- 
ness L = 20 nm has kp « 0.29A/a. For Py at 300 K, 
this gives the estimate Tims ~ 2.1 nm. If the thickness 
is reduced to 10 nm, then kp w OMOA/a and the RMS 
orbital radius increases to Tj-ms ~ 4.0 nm. Even though 
these are rather small, the distributions p{r) are rather 
wide and therefore at times one can expect even larger 
vortex gyrotropic oscillations. 

Finally we note that the thermal distribution of the 
vortex rotational velocity is connected to the radial dis- 
tribution p(r), because the Hamilton equations (j67|l im- 
ply 

V-^gxX, cog^-'^z. (77) 

Thus, we can transform magnitudes with V = ojgt- Then 
the RMS rotational velocity is 



K-ms =\^^g\ frms = ^ , (78) 

which varies proportional to \fkp jL. This is connected 
to a Boltzmann distribution for the probability /(F) dV 
of vortex speed Y in some interval of width dV ^ where 

/(y) ^ ^Jy^A ^ ^^^Ve-hP^oV\ (79) 

This involves a gyrotropic effective mass togj 

_G2^(2^)2 R 
"^«=fc;" 0:878^' 

determined both by the vortex force constant and by the 
disk thickness contained in the definition of G. For small 
aspect ratio, however, the thickness cancels and this mass 
is proportional to the disk radius alone. At i? = 100 
nm, the mass is about 1.2 x 10~^^ kg, independent of 
the material. Although /(V) has a mathematical form 
identical to that for p(r), it leads to another interesting 
interpretation of the vortex dynamics in equilibrium. 
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